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5-H . A numerical model is presented to describe both the transient and steady-state dynamics 

of viscous threads falling onto a plane. The steady-state coiling frequency is calculated as 
a function of fall height H . In the case of weak gravity, Vl oc and Vl (x H are obtained 
for lower and higher fall heights respectively. When the effect of gravity is significant, the 
. relation oc is observed. These results agree with the scaling laws previously predicted. 

O . The critical Reynolds number for coil- uncoil transition is discussed. When the gravity is weak. 



the transition occurs with hysteresis effects. If the plane moves horizontally at a constant 
speed, a variety of meandering oscillation modes can be observed experimentally. The present 
model also can describe this phenomenon. The numerically obtained state diagram for the 
meandering modes qualitatively agrees with the results of experiment. 
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0> ■ 1. Introduction 

T^lj- ■ The instability of viscous fluid generally occurs when the Reynolds number exceeds a 

CN ' critical value, because, in the low-Reynolds-number regime, the eigenmodes with short wave- 

00 ! length are impeded. However, when the fluid surface can move freely and deform greatly, this 
is not the case. For example, an axisymmetric jet emitted from a nozzle onto a horizontal 

:H ■ plane shows buckling instability below a critical Reynolds number. The "rope" of fluid loops 

;h ' near the plane and forms a coil owing to the buckling. This phenomenon is thus called "fluid 



rope coiling" and has been studied for several decades in the laboratory. Barnes and Wood- 
cock first performed an experimental study and observed that the coiling frequency increases 
proportionally to the fall height.^' ^ A more comprehensive investigation was conducted by 
Cruickshank and Munson,^ who found the existence of the critical Reynolds number above 
which bucking does not occur. They also found that the coiling frequency is not a mono- 
tonically increasing function of fall height but decreases for low fall height. The theoretical 
treatment of fluids with freely moving surfaces is not easy, however, some researchers have 
succeeded in determining the critical fall height and frequency at the onset of coiling using 
linear stability analysis with some simplifications and assumptions.^'^ 

* Multiple authors and affiliations correspond using arabic numerals each other. 
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Phenomcnologically, one can understand that steady coiling proceeds with a mechanical 
balance between the driving force of a steady flow and the internal viscous stress in the 
buckling portion of the rope.^ assumed that, in the high frequency limit, the inertial force Fj 
is equal to the viscous force Fy in magnitude. Representing the flow rate as Q, the radius 
of the rope in the coil as a, and the kinetic viscosity as this mechanical balance yields a 
scaling law for the coiling frequency, 

which is called "inertial coiling". Considering the gravitational force Fq and the fluid injection 
force Fp^ asserted that there are two more scaling laws. "Gravitational coiling" takes place 
under the condition Fq ~ Fy, and "Viscous coiling" when Fp ~ Fy - Each mechanical balance 
yields 

fi.oc(g)"\ (2) 

Sivoc^, (3) 

where g is the gravitational acceleration and H is the fall height. The existence of these three 
distinct coiling regimes were confirmed experimentally.^'^ 

Ribe also derived the differential equations for a very thin rotating rope to predict the 
steady-state coiling frequency as a function of fall height, and demonstrated the three different 
coiling regimes. His analysis revealed that the steady solution is multivalued so that there 
might a discontinuity in the selected coiling frequency.^' 

The theories of the buckling and coiling instability have been restricted to the description 
of the steady-state coiling frequency or the onset of buckling with infinitesimal amplitudes. 
Furthermore, a recent experiment revealed that fluid rope falling onto a moving belt shows a 
rich variety of "meandering" patterns, ^^"-"^^ but a theory that successfully explains the state 
diagram for various patterns has not yet been proposed. Hence a numerical model that can 
describe the entire dynamics of viscous fluid rope is still needed. 

In this study, we use our recently proposed numerical model^^ to understand the transient 
dynamics of the coil-uncoil transition and steady-state coiling frequency. We also extend the 
model to describe the meandering instability of fluid rope on a moving surface. 

This paper is organized as follows. In section 2, partial differential equations that describe 
the dynamics of coiling are derived. Coil-uncoil transition and its hystereric effects are dis- 
cussed in section 3. The coiling frequency as a function of fall height is discussed in section 
4. In section 5, the present model is applied to the problem of the meandering instability of 
fluid rope falling onto a moving surface. In section 6, we summarize our results. 
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2. Model 

We present a simple numerical model for a thread of viscous fluid falling onto a plane. 
Analogous equations were used by Chiu-webstar and Lister to illustrate the steady-state vis- 
cous catenary dragged by a horizontally moving belt.^^ 

The present model is restricted to the case of slight deformation: viscous stiffness for 
bending and twisting are neglected. We assume uniform flow within the rope and do not 
include the surface tension effect. Consider a viscous fluid ejected from an orifice at a sufficient 
height. We require that the upward growth speed of the coil be balanced by the downward 
slumping speed at the top of the coil, namely, the top of coiling portion does not move. Under 
this assumption, we fix the origin of the reference frame at the point where the fluid rope 
begins to coil (see Fig. 1) and only consider the flow at z > 0. 

Because of the assumption of slight deformation, the unit tangential vector of the rope 
t = {tx,ty,tz) should be almost parallel to the z axis. We thus replace the derivative with 
respect to t by that of z. Let S be the cross-sectional area parallel to the xy plane, and w the 
axial flow velocity. The conservation of the volume flux is written as 

Let q = {qx,Qy) be the center of mass at a given height, and u = {ux,Uy) be the velocity; 
these obey 

^ + ^^J<Z. = u., (5) 

where i = x and y. We denote the axial stress acting on the cross section as at, which would 
be the driving force of the oscillation, and the viscous shear stress as Uj, which would be the 
resistance force to bending. The equation of motion for w is 
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where p is the density of fluid. The ith component of the stress acting on a cross section is 
(Tj + tiGt- Therefore, we obtain the equation of motion for Ui as 

^ + "'l)"< = ^l'^ + « 

The shear stress can be written as (Xj = 'qd{ui + wti)/dz, with 77 as the viscosity. We can 
easily find that the axial stress at does not explicitly include the fluid pressure p, as follows. 
Let r be the radial coordinate for the local cross section, and Ur the flow velocity in the 
direction r; the conservation of volume flux gives 2{duj./dr) = —dw/dz. The radial stress can 
be written as Gr = —p + 2ri{du, /dr) = —p — dw/dz. The radial stress must vanish at the free 
surface, thus we obtain p = ~r]{dw/dz). The axial stress, therefore, can be written as^^ 

dw dw 

Using the expressions of at and ai, the dimensionless forms of Eqs. (6) and (7) become 

d d\ 3d f^dw\ 1 



dt dz J ^ SRe dz \ \ dz ^ dz dz 

(10) 

where Re = d\win\/v is the Reynolds number and Pr = wf^/gd is the Proude num- 
ber. Note that the vector t is determincxl from th(^ dcniAatiA^c of q with respect to z as 



n = {dq:c/dz,dqy/dz, 1) /C with C = ^{dq^c/dz)^ + {dqy/dzf + 1. 

Next, we discuss the boundary condition of the present model. At the neighborhood of the 
oriflce, we neglect the relaxation of Poiseuille flow to plug flow. Thus, at the injection point 
z = H/d, 

qi{H/d)=0, Ui{H/d) = 0, S{H/d) = ^, w{H/d) = -1. (11) 

At z = 0, we assume the free-end boundary condition 

g^(0) = 0. u^(0)=0, (12) 

where the prime indicates a derivative with respect to z. To determine the boundary value 
of the axial flow velocity w{0), we utilize a phenomenological parameter, the "energy loss 
coefficient (= a)" proposed by Cruickshank.^^ They postulated that the rope starts to buckle 
at the height z = ( where the viscous stress changes its sign. Considering energy, momentum 
and volume flux conservation across the buckling region < z < (, they derived w{0)/w{Q = 
(1— Q;)/(l-|-a). The value of a is experimentally determined as 0.76, which is fairly independent 
of viscosity, flow rate and oriflce diameter. Therefore, we require the following time dependent 
boundary condition: 

w{0) = PwiO, (5 = 0.14. (13) 
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Fig. 2. Trajectories of the bottom of the model rope starting at i = 0. (a) Trajectory when the 
circular coiling motion is stable in the steady state (Re=1.0, Fr=1.0, H/d=10.0). (b) Trajectory 
that converges to stable axial flow. Parameters are the same as in (a) but Reynolds number is 
slightly larger (Re=3.0). 




Fig. 3. Coiling regime in the Re-Fr plane for increasing (upward triangles) and decreasing (downward 
triangles) Reynolds number. The inset shows the radius of coiling for Fr — as a function of 
Re. 



Note that, because of Eq. (8), wlQ is the maximum axial velocity that can be determined by 
solving Eq. (9). How the present model is affected by the value of /? is discussed in the next 
section. 

We numerically solve the partial differential equations (4), (5), (9) and (10) using the Euler 
scheme. The space interval Az = 0.1 is fixed and the time step is set as At = 0.25ReAz^. 
Control parameters are Reynolds number Re, Froude number Fr, and fall height H/d. The 
simulation starts with the initial conditions w{z) = —1, s{z) = it/4:, qi{z) = and Ui{z) = 0. 
A small roughness with an amplitude of 0.01 is given to qi{z) as an initial shape at t = 0. 
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3. Coil-uncoil transitions 

It is empirically known that fluid threads do not buckle if the Reynolds number of the 
flow is sufficiently large. Experimentally, the critical value of Reynolds number for coil-uncoil 
transition ranges from 0.7 to 2.0 and in average 1.2, regardless of the value of Fr.'^ The present 
model correctly reproduces this well-known behavior. In Fig. 2(a), we plot the trajectory of 
q in the plane 2 = for Re = Fr = 1 and H/d = 10. The initial roughness of q{z) gradually 
increases and the trajectory converges to a circle with a radius of order unity. However, in the 
case of Re larger than the critical value (Re = 3), the amplitude of oscillation monotonically 
decreases, as shown in Fig. 2(b), and the axial stagnation flow becomes stable. 

We investigate critical Reynolds number Re*, and find that the transition occurs with 
a hysteresis effect. In the simulation, we change the value of Re sufficiently slowly (slower 
than 0.1% change per cycle of motion) after the motion of fluid reaches a steady state. A 
small perturbation of the amplitude, lO^^d, is continuously added for qi throughout the time 
evolution. In this way, we obtain the coiling radius as a function of increasing or decreasing Re, 
as shown in the inset of Fig. 3. For increasing Re, the radius sharply falls at Re = 1.6 (= Re^) 
and for decreasing Re, the radius rises at Re = 1.0 (= Re^*). The axial stagnation flow is 
absolutely stable at Re > Re^ and is absolutely unstable at Re < Re^*. In the bistable state 
(Re^ < Re < Re^), the fluid thread starts to coil if the amplitude of perturbation exceeds 
unity in order. The difference Re^ — Re^ reaches the maximum at Fr = 10^ and disappears 
for Fr < 10. Although the hysteresis effect in the coil-uncoil transition has not been observed 
experimentally, such behavior can occur owing to the inertial effect (centrifugal force keeps 
the fluid thread rotating). 

Next we show that the phenomenological value /? = 0.14 also can be justified by comparing 
the critical Reynolds number. The experiment performed by Cruickshank^ yielded the critical 
value with increasing flow rate, and it corresponds to Re^ in the present model. As shown 
in Fig 4(a), Re^ strongly depends on /3, however (3 = 0.14 gives Re^ = 1.38 which is in 
accord with the experimental value Re^ = 1.2. In contrast to the critical Reynolds number, 
the dependence of coiling frequency on /3 is relatively slight as shown in Fig 4(b). Thus the 
discussions in the next section are not strongly influenced by the choice of (5. 

4. Frequency of steady state coiling 

4.-1 Numerical results 

In this subsection, we show coiling frequency as a function of fall height. The experiments 
reported by Habibi et al.^ are performed under Re ~ 10~^, which requires very short time 
step for computation. The numerical results shown below are obtained under Re around unity, 
thus we mention that the comparisons with the experiments are qualitative. 

Figure 5(a) shows the steady-state coiling frequency as a function of fall height under 
weak gravity conditions (Pr » 1). The frequency decreases with slope —1 for H/d < 4, and 
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Fig. 4. (a) Critical Reynolds number RejJj versus f3. (b) Steady-state coiling frequency as a function 
of j3 for Re = 0.4. Other parameters are Fr = 10* and H/d = 10 in both plots. 




Fig. 5. Coiling frequency as a function of fall height. Insets is the experimental result by Maleki et 
al.^ (a) Coiling under weak gravity with Re= 3, and Fr= 100 (circles), Fr=200 (triangles) and 
Fr=400 (squares), (b) Coiling under strong gravity with Re=: 0.1 and Fr= 0.01 (circles), Fr=0.02 
(triangles) and Fr=0.05 (squares). 



becomes almost constant around H/d = 10. For higher fall height, the influence of gravity 
becomes significant, and the frequency shows a linear increase. The inset shows the same plot 
as obtained experimentally, which is in good agreement with the present model. 

Next we show the coiling frequency under the strong gravity condition (Fr ^ 1) in Fig. 
5(b). Here, we choose a characteristic time scale \J dj g. In this case, the viscous coiling regime 
is negligible, and the frequency strongly depends on Fr. The frequency-height curve has the 
slope of 2.0 while the experiment shows the slope of 2.5. A clear discrepancy appears for lower 
fall height, at which the experiment shows a marked discontinuous jump. The present model 
does not reproduce this behavior. 
scaling laws 

To compare frequency-height relationships with the the scaling laws [Eqs.(l-3)], we must 
estimate the rope radius a, which depends on other flow conditions {Q, g, v and }i\ Let C, be 
the height at which the flow velocity w reaches a maximum. Under the action of gravity, the 
internal stress of the fluid rope changes from tensile to compressive at this point, and buckling 
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must occur at z < Therefore, we reasonably assume a to be the radius at z = and divide 
the fluid rope into the "tail region" {z > (") and the "buckling region" (z < ("). 

Radius a is governed by the gravity-induced thinning which is determined by the steady- 
state form of Eq. (6): 

ww' = Svw ^ — ^ — g, (14) 

To solve this one-dimensional problem, we consider only the flow in the tail region and set 
z = ( as the origin of the z axis. The boundary conditions are w{H) = —1 and w'{0) = 0. 
First, we consider the case of weak gravity. For g = 0, Eq. (14) has an obvious solution, 
w{z) = —Win- Thus, we seek a solution with the form 

w{z) = -Win + g4'{z). (15) 

Substitution of Eq. (15) into Eq. (14) yields 

3u(f)"{z) + Wi^(l)'{z) -1 = 0. (16) 

Solving this equation with the boundary condition i^'(O) = 4>{H) = 0, we obtain the pertur- 
bation solution 

^ = -1 + 1 f ^ + e-/-^ - e-""^) , (17) 
Win Ft V d J 

where we set 3/Re = 1 for simplicity. The conservation of volume flux gives the radius a as 



^-^l + ^f^.e-Z-lH"" (18) 



^ w{Q) ' Ft V d 

Assuming H/d ^ Fr, Eq. (18) implies that a ^ d. Under this condition, the scaling law 
of viscous coiling may take place, therefore, we obtain fly oc Qd~^H~^. As the fall height 
increases, the gravity becomes significant for the tail region. For H/d !=a Fr » 1, Eq. (18) 
can be approximately written as a = d{l + gH /w\^)~^f'^ . Substituting this relation into the 
scaling law (2), we obtain 



This indicates that the frequency increase proportional to iJ, which are in accordance with 
the numerical result shown in Fig. 5(a). Thus, the frequency that changes from decrease to 
increase corresponds to the transition from viscous to gravitational coiling. 

Next we consider that both Re and Fr are much smaller than unity. In this case, the 
inertia term in Eq.(14) is negligible, thus w{z) is determined by the balance of gravitational 
and viscous forces. For the higher fall height, the flow is strongly stretched (a ^ 1). In this 
limit, the force balance implies a ~ {Qu/gy/^H"^.^ In the buckling region, both gravitational 
and inertial forces are important. Eqs. (1) and (2) with the strong stretching condition yield 

/ „5 / 5 \l/3 
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Fig. 6. Frequency versus 1/Fr for Re — 0.1 and _ff = 30 (circles), H — 15 (triangles) and H ^ 4 
(squares). Each solid line has the same slope of 1.25. 
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Fig. 7. (a) Numerically calculated value of C for Re=: 3.0 and Fr= 100. (b,c) Viscous-gravitational 
transition height H* as a function of 1/Fr and \/Re, respectively. 



Because the H has similar exponents in these two relations, the gravitational and inertial 
coiling may behave similarly with changing fall height. Experimentally, oc was reported 
by Cruickshank and Munson.'^ Maleki et al. also reported Q. oc H"^'^.^ To judge whether 
these results obey the scaling law of gravitational or inertial coiling, we examine the gravity 
dependence of frequency. In Fig 6, we show frequency as a function of 1/Fr for Re = 0.1 and 
Fr > 10. Each plot lies on the slope 5/4. We therefore conclude that the frequency shown in 
Fig. 5(a) corresponds to gravitational coiling. 

The present model does not reproduce the scaling law of inertial coiling, or the discontin- 
uous jump of frequency observed in the gravitational-inertial transitional regime. A possible 
root of this discrepancy is the boundary condition at the bottom of the fluid rope. In the 
experiment for the high-frequency coiling regime, the rope rapidly piles up to form a column 
and it collapses when the height exceeds a critical value. ^ In this case, the top of the coiling 
portion is not stationary, and relation (13) may fail. Therefore, the boundary condition at 
z = must be improved to render the present model applicable to the inertial coiling regime. 
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4.3 Viscous-gravitational transition 

Within the viscous coihng regime, the coihng frequency decreases with faU hight, while 
the frequency increases in the gravitational coiling regime. We discuss the critical fall height 
H* at which the frequency is minimized, and which corresponds to the transition between the 
viscous and gravitational regime. In figure 5(a), the transition occurs at H* ~ 6d. However, 
H* depends on viscosity and gravitational acceleration. In Fig. 7(a), we plot ( as a function 
of H/d for the weak gravity condition (Re = 3 and Fr = 10^). For shorter fall height, the 
maximum w appears at the point of injection (plots lie on the line C = H/d). Comparing C, 
with the frequency Qd/wi^ in Fig. 5(a), we find that the frequency decreases with H only 
when the relation ( = H/ d is satisfied. Therefore, viscous coiling appears only when the rope 
is wholly compressed. Because this feature can be seen for a wide range of Re and Fr as long as 
the viscous coiling regime exists, we claim that the H* is identical to the maximum of Q in Fig. 
7(a). The maximum of Q can be realized as a relaxation length of which the effect of boundary 
condition at z = can travel through a rope. We can uniquely construct a dimension of length 
using g and u, as g~^l^v'^l'^ . Therefore, 



We show iJ* as a function of l/Pr (Fig. 7(b)), and 1/Re (Fig. 7(c)). Each plot respectively 
lies on the slope -0.39 and 0.68 which agrees with the Eq. (21). 

5. Coiling and meandering of dragged fluid rope 

In this section, we apply the present model to a fluid thread falling onto a horizontally 
translating surface with speed U<^y^ experimentally found that the thread of fluid deposited 
on the moving surface shows a rich variety of meandering patterns due to the buckling of 
the rope. If the speed of the surface is sufficiently high, the falling thread is strongly dragged 
to form a viscous catenary. In this case, a straight line is drawn on the surface. As [/q 
decreases, the fiuid thread starts to oscillate and sinusoidal or other period-doubling curves 
appear. Finally, the oscillation reverts to the rotational one of the ordinary fiuid rope coiling 
onto a stationary plate, thereby a cycloid pattern appears on the moving surface. 

Although the present model does not include a bottom surface, we can introduce a hori- 
zontal drag by adding a control parameter to the boundary condition of Eq. (12) as 



Here, 7 is the shear rate emerging as a result of the surface motion. We mention that the 
shear rate may vary during the motion, however, in the present model, we omit that effect 
and regard 7 is positive constant. Hereafter, we fix Re = 0.1 and Pr = 3.2 x 10^^ (strong 
gravity condition). 

Changing the value of 7, we found that the thread exhibits three distinct states: translated 
coiling (TC), meandering (M) and catenary states. Figure 8 is the trajectory of the bottom of 



iJ* = max(C) oc 5-^/3^,2/3 



(21) 



<(0)=7- 



(22) 
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Fig. 8. Steady state trajectories of the bottom of the dragged fluid rope for H = 10. (a) Trajectory of 
translated coil state (s = 4). (b) Trajectories of the meandeing state (s = 10). (c) A intermediate 
state between translated coiling and meandering state (s=6). 




Fig. 9. State diagram in H-^ plane. Upward (downward) pointing triangles indicate the threshold 
for increasing (decreasing) shear rate. 



the thread under steady-state oscillation. As 7 increases from zero, circular coiling is gradually 
strained [Fig. 8(a)], and subsequently, that the trajectory changes to a figure-eight shape 
pattern [Fig. 8(b)]. If the shear rate exceeds a critical value 7c, the oscillation diminishes in 
amplitude and the thread forms a stationary catenary (not shown). 

Fig. 9 is a diagram of the various oscillation modes in the H-^ plane. The upward (down- 
ward) triangles indicate the boundaries of the modes for increasing (decreasing) shear rate. 
A very weak hysteresis effect can be seen only for lower fall height. The experimentally ob- 
tained state diagram is given by^^ and more precisely by.^^ Here, we note that, because the 
relationship between Uq and 7 is unknown, the present model cannot reproduce the patterns 
of a thread laid down on a moving surface, while the experiments categorize the modes of 
oscillation on the basis of the patterns. In this study, we compare Fig. 9 with the state diagram 
reported by^^ by assuming that 7 is proportional to C/q, and claim that the trajectories in 
Figs. 8(a) and (b) respectively correspond to TC and M states. 

If the fall height is lower than /7 ~ 8, the transition from TC to M state occurs with 
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Fig. 10. Amplitude in y direction and the ratio of frequency U,x/^y as a function of increasing 7. 
Fall height H/d^ 11. 

discontinuous jump in the amplitude in the y direction. However, for H > 8, the discontinuity 
vanishes and a crossover of TC and M states starts to appear. In this region, complex oscilla- 
tions, tje path of which is not closed, are observed, as shown in Fig. 8(c). It is notable that, 
in the experiment, M and TC states rarely arise for H > 8.7 ± 0.5 and complex oscillations 
named "stretched coiling" or the "W state" are observed for higher fall height. 

In Fig. 10, the solid line shows the oscillation amplitude in the y direction with increasing 
7. The amplitude gently increases with 7 in the TC region (7 < 6), and decreases in the 
M region (8 < 7 < 7c). Immediately below 7 = 7c, the amplitude rapidly falls to zero, as 
experimentally observed. In the same figure, we also plot the ratio of x to y components of 
the steady-state frequency Q^/^y During the transition from TC to M, This ratio jumps 
from 1/1 to 2/1, which agrees with the result of the experiment. We mention that, in the 
intermediate region (7 ~ 8), the ratio 0,x/^y slightly deviates from 1/1. 

Experimentally, it is known that increases linearly with 7 for 7 < 7c, and the frequency 
at the bifurcation point is in excellent agreement with that of the ordinary coiling state on 
the stationary surface. However, the present model does not reproduce this tendency; instead, 
the frequency 0,y is a monotonically decreasing function of 7. 

6. Summary 

In this paper, we have presented a simple numerical model for the motion of a viscous 
thread falling onto a plane. We found that the critical Reynolds number Re* for coil-uncoil 
transition strongly depends on the boundary condition [Eq. (13)] at the bottom of the thread. 
However, the phenomenological parameter /?, which was experimentally determined, gives 
good estimate of Re* . 

The coiling frequency and the scaling laws have been reviewed by classifying the parameter 
space into areas with weak and strong gravity. In the former, the present model is in accord 
with both the scaling laws and the results of experiments. However, in the latter case, the 
inertial regime cannot be reproduced. This may be because, in high-frequency coiling, the 
boundary condition for z = is no longer valid. We must overcome this difficulty in our 
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future work. 

The present model can describe the meandering instability of a viscous thread falling onto 
a moving surface simply by changing the boundary condition at the bottom of the thread. The 
obtained diagrams of meandering motion and amplitude are in qualitative agreement with the 
result of experiments. 
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